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ABSTRACT 

We explore a ballistic orbit model to infer the gravitational drag force on an accreting 
point mass M, such as a black hole, moving at a hypersonic velocity vo through a 
gaseous environment of density po. The streamlines blend in the flow past the body 
and transfer momentum to it. The total drag force acting on the body, including the 
nonlinear contribution of those streamlines with small impact parameter that bend 
significantly and pass through a shock, can be calculated by imposing conservation 
of momentum. In this fully analytic approach, the ambiguity in the definition of the 
lower cut-off distance r m ; n in calculations of the effect of dynamical friction is removed. 
It turns out that r m j n = yfeGM /2vq. Using spherical surfaces of control of different 
sizes, we carry out a successful comparison between the predicted drag force and the 
one obtained from a high resolution, axisymmetric, isothermal flow simulation. We 
demonstrate that ballistic models are reasonably successful in accounting for both the 
accretion rate and the gravitational drag. 
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1 INTRODUCTION 

A body moving in a background medium loses momentum 
due to its gravitational interaction with its own gravitation- 
ally induced wake. This process is often referred to as dy- 
namical friction. Chandrasekhar (1943) estimated the dy- 
namical friction on a massive particle passing through a ho- 
mogeneous and isotropic background of light stars. In the 
case where the perturber moves in a gaseous medium, the 
gravitational drag is traditionally inferred as the gravita- 
tional attraction between the perturber and its own wake. 
In this approach, the density structure of the wake is de- 
rived in linear perturbation theory by assuming that the 
body produces a small perturbation in the ambient medium 
(Dokuchaev 1964; Ruderman & Spiegel 1971; Just & Kegel 
1990; Ostriker 1999; Kim & Kim 2007; Sanchez-Salcedo 
2009; Namouni 2010). For a perturber moving on a rectilin- 
ear orbit at constant velocity, the steady-state linear theory 
predicts that the drag force vanishes for subsonic perturbers, 
while it becomes similar to the collisionless drag force for su- 
personic bodies. Ostriker (1999) considered the linear-theory 
drag as a time-dependent rather than steady state problem 
and arrived at the following formula for the gravitational 
drag force, 
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The perturber of mass M, which moves at velocity vo and 
Mach number M in a rectilinear orbit through a homo- 
geneous medium with density po and sound speed Co, is 
assumed to be formed at t = 0. The minimum radius 
r m ij, is the typical size of the perturber. This formula has 
enjoyed widespread theoretical application (Narayan 2000; 
Escala et al. 2004; Kim 2007; Conroy & Ostriker 2008; 
Villaver & Livio 2009; Tanaka & Haiman 2009; Nejad- 
Asghar 2010; Chavarria et al. 2010). Because of the linear- 
theory assumption, the above equation is properly valid 
only at r > -Rbh where Rbh is the Bondi-Hoyle radius 
(Rbh = GM/[4(1 + M 2 )]). Therefore, Equation (1) is 
strictly valid for extended perturbers with a softening radius 
much larger than the Bondi-Hoyle radius. In fact, for ex- 
tended perturbers, Sanchez-Salcedo & Brandenburg (1999) 
found good agreement between the gravitational drag in 
full hydrodynamical simulations and Ostriker's formula. In 
particular, for Plummer perturbers with softening radius r s 
much larger than Rbh, they found that r m i n = 2.25r a . An 
extension of Ostriker's formula for extended bodies orbiting 
in a stratified gaseous sphere was given in Sanchez-Salcedo 
& Brandenburg (2001). 
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In the case of point-like perturbers, like massive black 
holes, it is reasonable to assume that r m i n should be of the 
order of a few _Rbh, but a nonlinear analysis is required to 
fix the uncertainty in the definition of r m i n . In adiabatic sim- 
ulations of axisymmetric accretion flows past a gravitating 
absorbing object, Shima et al. (1985) computed the drag by 
considering two contributions: the aerodynamic force, which 
is due to the accretion of momentum over the body surface, 
and the gravitational force on the perturber by its own wake. 
They found that the numerical results were consistent with 
the estimates in linear theory. 

The problem of the gravitational drag on a point-mass 
particle has revived new interest to estimate the timescale of 
the orbital decay of massive black hole binary in the centre 
of galaxies. Escala et al. (2004) simulated the orbital decay 
of a single black hole moving initially on a circular orbit in 
an isothermal gaseous sphere. They found that the gravita- 
tional drag is less peaked at 1 < A4 < 2 than predicted by 
Ostrikcr's formula with ln-uo£/?" m m = 3.1. Tanaka & Haiman 
(2009) combined the prescriptions of Ostriker (1999) and Es- 
cala et al. (2004) into a formula that is used as a prescription 
of the gaseous drag on black holes in numerical simulations. 
In order to isolate the physical reason of the failure of Os- 
triker's formula, Kim & Kim (2009) and Kim (2010) car- 
ried out axisymmetrical simulations of a massive body in 
rectilinear orbit with different values of the strength of the 
gravitational perturbation due to the body as measured by 
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where r s is the softening radius of the Plummer perturber. 
They find that the functional form of the gravitational drag 
is not so peaked as the linear theory predicts and conclude 
that the discrepancy between the numerical and Ostriker re- 
sults are most likely due to the nonlinear effect. It is impor- 
tant to note that in the simulations of Escala et al. (2004), 
Kim & Kim (2009) and Kim (2010), the perturber simply 
provides a smooth gravitational potential and does not hold 
any absorbing surface. Without any absorbing inner bound- 
ary condition, a hydrostatic envelope with front-back sym- 
metry is formed near the perturber. Because of the front- 
back symmetry, this large envelope provides a negligible con- 
tribution to the gravitational drag force. 

However, it is well-known that accretion is a crucial in- 
gredient in point-like objects, such as black holes or stars. 
Different boundary conditions in the high density region of 
the wake are expected to change the gas dynamics near the 
perturber and the strength of the gravitational drag (e.g., 
Fryxell et al. 1987; Naiman et al. 2011). For example, Ruffert 
(1996) simulated a 3D quasi-isothermal flow with an absorb- 
ing boundary surrounding the point-like object. Moeckel &: 
Throop (2009) carried out similar simulations, but then also 
included an accretion disk orbiting the point source. 

The aim of this paper is to describe the contribution 
of the nonlinear inner wake to the gravitational drag on 
hypersonic perturbers by using the ballistic orbit theory 
(Bondi & Hoyle 1944; Lyttleton 1972; Bisnovatyi-Kogan et 
al. 1979). Whereas this theory (the so-called model of line- 
accretion) has been extensively used as a powerful frame- 
work to describe the gravitational interaction between a 
moving massive body and the surrounding gaseous medium 
in the context of supersonic Bondi-Hoyle- Lyttleton accretion 



(e.g., Koide et al. 1991; Edgar 2004), it has been tradition- 
ally ignored as a tool to quantify the gravitational drag. In 
fact, all analytical studies about the gravitational drag in 
gaseous media have been based on the linear perturbation 
theory, following on the analysis of Dokuchaev (1964), Rud- 
erman & Spiegel (1971) and Rephaeli & Salpeter (1980). In 
this paper we develop the ballistic orbit theory to provide an- 
alytical expressions of, not only the mass accretion rate, but 
also the nonlinear drag force on a hypersonic compact body. 
These estimates will be compared with numerical results of 
an axisymmetric isothermal hydrodynamical simulation. 



2 THE FREE-STREAMING FLOW SOLUTION 

Let us consider the axisymmetric flow generated by a point 
mass M which moves hypersonically at a constant veloc- 
ity vq inside a homogeneous gaseous environment (see also 
Bisnovatyi-Kogan et al. 1979) . Fig. 1 shows a schematic di- 
agram illustrating the trajectory of a fluid parcel in a frame 
of reference at rest with respect to the point mass. 

Because vo (the upstream environmental velocity, see 
Fig. 1) is hypersonic, we neglect the pressure force and con- 
sider the ballistic trajectory of the fluid parcels in the gravi- 
tational potential of the point mass. As they have a positive 
E — v%/2 energy (per unit mass), the trajectories of the 
fluid parcels are hyperbolae of the form : 



£o (1 + cos (9) sin (9 ' 
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where £ is the impact parameter of the fluid parcel (see Fig. 
1) and 
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with G being the gravitational constant. For deriving Eq. 
(3) one has to consider a generic hyperbolic trajectory of the 
form r = p/[l + e cos(0 — 0o)], and then impose the upstream 
boundary condition and the conserved angular momentum 
£«o to determine the constants p, e and 6o- 

From Eq. (3), we see that the streamline intercepts the 
symmetry axis (i.e. — 0, see Fig. 1) at a position 



xo 



2£o ' 
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downstream from the perturber. The material will therefore 
pile up in a narrow, downstream wake surrounding the sym- 
metry axis, forming a dense column of gas. 

From the equation for the streamlines (Eq. 3) one can 
calculate the velocity components of the free-streaming flow 
along the x and j/-axes : 



v x — ^ (£ + Co sin 6) ; v y 



(l + cos<9) 



(6) 



Assuming that the environment has a homogeneous density 
po far upstream from the source, it is possible to obtain the 
density p(x,y) of the free-streaming flow as a function of 
position: 



p _ e e 

with 



(7) 



Gravitational drag 3 







— ■ - 






M x 



Figure 1. Schematic diagram showing the trajectory of an en- 
vironmental fluid parcel in hypersonic motion with respect to a 
point mass M. The initial velocity vn of a parcel with impact pa- 
rameter £ is parallel to the rr-axis, and its trajectory is the 7(8) 
curve. The problem has cylindrical symmetry, with x being the 
symmetry axis and y the cylindrical radius. 



\/x 2 + y 2 ; £ = i [y + Vv 2 + 4£o(r + x)] ■ 



(8) 



Note that x is the distance along the symmetry axis and 
y the cylindrical radius. It is simple to see that p > pn. 
Equation (7) is not valid in the shocked column of gas near 
the positive £-axis. The density enhancement in this ap- 
proach is different from that derived in linear theory which 
predicts zero-enhancement outside the Mach cone (e.g., Os- 
triker 1999). Koide et al. (1991) found a good accordance 
between the analytical solutions (6)-(8) and the numerical 
solution even for Mach numbers as low as 1.4. 

From Eq. (6) we see that at the point in which the 
streamlines intercept the symmetry axis (i.e., for 8 = 0), 
the flow velocity has an ^-component v x (0) = vo (identical 
to the far upstream flow velocity and independent of the 
impact parameter £ of the flow parcel) and a {/-component 
v »(0) = ~2t>o£o/£- This latter component of the velocity 
will be thermalized in a shock surrounding the downstream 
wake. We will assume that the post-shock thermal energy is 
radiated away instantaneously. 

Now, the kinetic+potential energy per unit mass of the 
flow at x — y — oo is En = Vo/2. When the flow hits the 
symmetry axis (at = 0), the energy associated with the y- 
velocity is thermalized, so that the kinetic+potential energy 
is reduced to a value 



Et — En 
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From Eq. (9) it is clear that E t < (i.e., the post-shock 
material is gravitationally bound) if the condition 



(10) 



is met. Therefore all of the material arriving with impact 
parameters < 2£o will eventually be accreted onto the body. 
A streamline with impact parameter £i = 2£o (with £o given 
by Eq. 4) crosses the symmetry axis at a distance 

xi = 2£o , (11) 

downstream from the body (see Eq. 5 and Fig. 1). 

The material within the downstream wake will have a 
complex flow pattern. From Eq. (6) it is clear that the ma- 
terial enters the tail with a positive rc-velocity. The material 
with impact parameter £ < 2£o (which is gravitationally 
bound, see above) will therefore enter the wake flowing in 



the -l-x-direction, so that it will first flow away from the 
body, and eventually reverse and fall back onto the body. 
The distance x m from the body (along the a;-axis) at which 
the flow reverses can be obtained from the condition of zero 
velocity for a radial motion in the gravitational potential. 
This condition gives : 
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We see that when f — y £i, x m — y oo. Consequently, stream- 
lines with impact parameter close to and smaller than £i, 
will take a long time to be accreted. 

The material in the wake that remains gravitationally 
unbound when entering the wake (i.e., the material with 
impact parameters £ > 2£o, see above), will flow away from 
the body along the x-axis, reaching infinity with a velocity 



2£o 



1/2 
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3 THE MASS ACCRETION RATE AND THE 
DRAG FORCE 

3.1 The accretion rate 

From the solution of section 2, the accretion rate onto the 
point mass can be obtained (Hoyle & Lyttleton 1939; Bondi 
& Hoyle 1944). All of the material with impact parameters 
£ — 2£o (see eq. 10) will fall onto the body. Therefore, the 
mass accretion rate is: 

4ir(GM) 2 po 



M acc = 7t(2£ ) Po«o = 



(14) 



where we have used Eq. (4) for the second equality. 



3.2 The gravitational drag 

We calculate the gravitational drag on the perturber by com- 
puting the net ^-momentum per unit of time, going 
through a spherical control volume of radius R > 2£o cen- 
tred on the body. It is clear that the contribution to the drag 
by the gas within the sphere is equal to Il x . 

We consider the three streamlines shown in the 
schematic diagram of Fig. 2 : 

• a streamline with impact parameter £i = 2£o, which 
crosses the axis at xi = 2£o (see eq. 11). All of the material 
with £ < £i is accreted onto the body, 

• a streamline with impact parameter £2 = v / 2£o~R, which 
crosses the axis at a distance R downstream from the per- 
turber (where R is the radius of the control sphere, see 
above) , 

• a streamline with impact parameter £3, which tangen- 
tially touches the control sphere (at a point with polar angle 
03, see Fig. 2). 

In order to obtain £3 and 63 we first set r = R (i.e., a 
radius equal to the radius of the control sphere) in Eq. (3), 
and invert this equation to find : 



sin 9± 



ZoR±Zoy/W< 
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which gives the two values of 9 at which the streamline with 
impact parameter £ cuts the control sphere. The angle 8+ 
(obtained with the + sign of the right hand term of Eq. 15) 
corresponds to the point in which the streamline enters the 
control sphere, and 6- corresponds to the exit point. For the 
tangential streamline (with impact parameter £3, see Fig. 2) 
the entry and exit points coincide, so that the term within 
the square root of Eq. (15) has to be equal to zero. From 
this condition, we obtain the value of £3 : 



(16) 



(17) 



£3 = v/i?(2& + R) , 
and using Eq. (15) we then obtain 

. yfl(2£ + R) 

Sm0a= fl + ft ■ 

Now, the i-momentum entering the control sphere from 
the upstream region can be calculated as : 

tl x , in (R) = 2tv Po v v x (£, 6+) £,di , (18) 
Jo 

where £3 is given by Eq. (16), v x (£,d) by Eq. (6) and 9+ is 
obtained with the + sign of Eq. (15). This integral can be 
solved analytically to obtain : 



Hx,in{R) 
where 
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1 + W() 



l + w a + VI + 2wo 



,(20) 



with w = £o/R < 1/2. 

The i-momentum rate leaving the control domain has 
two terms: 

• the rate fl Xt b of i-momentum leaving through the 
boundary of the spherical domain, 

• the rate tl x , a of ^-momentum hitting the symmetry axis 
and exiting the domain through a narrow wake along the x- 
axis. 

The momentum rate leaving the sphere through the 
boundary of the control volume is given by : 

U x , b = 2-Kpov [ 3 v x (£, 9-) £d£ . (21) 
J £2 

This integral can be performed analytically to obtain : 
tl x ,b = 4n^oPoVo fb , (22) 
where 
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(23) 



The momentum rate exiting the control region through 
the wake is given by : 



n*,a = / 2 



vr dm , 



where 

dm = 2-K^poVodt; , 



(24) 



(25) 



and the velocity vr along the axis with which the material 
leaves the control domain is given by the kinetic+potential 
energy conservation condition 



^o 2 



GM 

Xo 



GM 
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(26) 



where xo is the distance along the a>axis at which the 
streamline intercepts the axis as given by Eq. (5). Here we 
have used that the ^-component of the velocity is vo, as de- 
rived in Eq. (6). Using Eqs. (25-26), the integral in Eq. (24) 
can be performed analytically to obtain : 



tlx, a = 47r£oPo«< 2 fa ■ 

where 



(27) 



_ 1-(2 W0 ) 3 / 2 
U ~ 2w + 
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2w + ^2^(1 + 2w ) 
1 + VI + 2w a 



(28) 



Finally, the net drag force of the gas on the body is 
obtained as : 

F d = ii x , in - ii x , b - ii x , a = 4 ^ (Gi f )2p0 R) , (29) 



where £0 is given by Eq. (4) and 

f(t0,R)= fin- fb~ fa, 



(30) 



with fi n , fb and f a given by Eqs. (20), (23) and (28), re- 
spectively. Note that Equation (29) includes the force on 
the body due to momentum accretion (sometimes called as 
aerodynamic force). 

It is straightforward to show that in the R 3> £0 limit 
(too <§C 1, sec Eqs. 20, 23 and 28), this function takes the 
form 



2R 

e 



(31) 



A comparison between the full (Eq. 30) and approximate 
(Eq. 31) forms of / is shown in Fig. 3. It is clear that for 
radii larger than ~ 5£o the two forms of / agree to better 
than ~ 10% . 

From Eq. (31) we see that the drag force diverges loga- 
rithmically for large values of R, as occurs in linear pertur- 
bation theory. Therefore, it is possible to match the solution 
found in linear theory (Eq. 1), which is valid at far enough 
distances from the body, with that found in the nonlinear 
analysis. This can be accomplished by replacing R for vot, 
where t = is the time at which the body is formed 1 . More- 
over, we see that the ambiguity in the definition of the mini- 
mum radius r m i n that appears in linear theory is removed in 
our framework. In fact, we find that r m i n = Ve£o/2, where 
the factor Ve comes from inserting the term —1/2 that ap- 
pears in the right-hand-side of Equation (31) in the argu- 
ment of the log. 



1 In a realistic situation, R will increase with time until it reaches 
the boundary of the cloud. 
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Figure 2. Schematic diagram showing the control sphere (of ra- 
dius R > xi = 2£o ) and the three streamlines used in the calcu- 
lation of the drag force (see section 3.2). 
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Figure 3. Exact (solid line, sec Eq. 30) and approximate (dashed 
line, see Eq. 31) forms of the / function, which gives the depen- 
dence of the drag force as a function of the radius R of the control 
volume. 



4 AN AXISYMMETRIC NUMERICAL 
SIMULATION 

We have computed an axisymmetric numerical simulation, 
solving the Euler equations for an isothermal flow in a uni- 
form, cylindrical computational grid. We have used the "flux 
vector splitting" algorithm of van Leer (1982), with the sec- 
ond order (time and space) implementation described by 
Raga et al. (2000). 

The simulation that we are presented can be compared 
with the work of Ruffert (1996) and Moeckel & Throop 
(2009) , who computed 3D simulations of basically the same 
physical situation. The main difference with this previous 
work is that our simulation is 2D (axisymmetric), and has 
~ 2 orders of magnitude higher resolution. 

The computational domain has an axial extent of 15£o 
(with £o being the gravitational radius given by Eq. 4) and 
a radial extent of 7.5£o, resolved with 9000 x 4500 (axial x 
radial) grid points. A point mass (influencing the flow only 
through its gravitational attraction) is placed in the middle 
of the axial extent of the domain. A spherical volume of 



radius 0.05£o (30 pixels) around the body is artificially kept 
at a low density at all times, so that the material entering 
this volume from the rest of the computational domain is 
effectively removed. We simulate in this way the accretion 
of gas onto the object. 

In the left boundary of the domain we impose an in- 
flow of density po and velocity vo , parallel to the symmetry 
axis. A reflection condition is applied on the symmetry axis, 
and a zero gradient condition is applied in the remaining 
two boundaries of the computational domain. In the initial 
condition, the domain is filled with a uniform flow (of ve- 
locity vo and density po) parallel to the symmetry axis. The 
isothermal sound speed is chosen to be Co = «o/5 (i.e., the 
flow entering the domain has a Mach number of 5). 

The results obtained after time-integrations of 10, 40 
and 70£o /«o are shown in Fig. 4. This figure is a zoom 
of an inner region of the computational domain, showing 
the highly time-dependent wake formed downstream of the 
body. 

We take the density and flow velocity time-frames ob- 
tained from the simulation, and compute the net mass Mace, 
and momentum fluxes through a control sphere of arbitrary 
radius R centred on the point mass. Thereby, only M acc 
computed with R = 2£o corresponds exactly to the mass 
accretion rate onto the body. We also compute the gravi- 
tational force exerted on the body by the material within 
the control volume. In Fig. 5, we show the mass flux M acc , 
the drag force Fd and the gravitational force F g computed 
with a control volume of radius R = 5£o- F g is inferred as 
the gravitational attraction between the body and the per- 
turbed medium. M acc and Fd show a peak at t « 5£o/«o, 
and have fluctuating values for t > 10£o/^o- The gravita- 
tional force F g initially grows as more material enters the 
wake behind the object, and also shows fluctuating values 
as a function of time. The fact that M acc , Fd and F g have 
strong fluctuations is not surprising given the strongly time- 
dependent structure of the flow (see Fig. 4). 

In order to carry out a comparison with the analytic 
model (see section 3.2), we have calculated the average val- 
ues and the dispersions of M acc , F d and F g in the interval 
10£o/vo < ^ < 66£o/^o, in which the fluctuations of these 
quantities appear to be statistically stationary (see Fig. 5). 
We then plot these time-averaged values as a function of the 
radius R of the control volume in Fig. 6. 

We have considered control volumes with radii 2£o < 
R < 7£o, the lower boundary being fixed by the derivation 
of the analytic model (in which it was assumed that R > 2£o, 
see section 3) and the upper boundary given by the approach 
to the outer edge of the computational grid. It is clear from 
Fig. 6 that the dispersions of the M acc and Fd values grow as 
a function of R (due to the fact that larger, more massive ed- 
dies are seen at larger distances downstream from the body, 
see Fig. 4), while the dispersion of F g (which is a quantity 
integrated over the volume of the control sphere) remains 
approximately constant. 

The analytic model predicts that M acc = 47r£oPofo (see 
Eq. 14) for all control spheres with R > 2£o- The top panel 
of Fig. 6 shows that the time-averaged values obtained from 
the numerical simulation closely reproduce this result. The 
central panel of Fig. 6 shows the drag force F d calculated 
using Eq. (29), which has values that differ from the results 
from the numerical simulations by less than ~ 15 % (though 
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Figure 4. Density (in units of po, with the colour scale given by the top right bar) and velocity field (white arrows) from the axisymmctric 
simulation described in section 4, obtained for integration times t = 10f;o/ w O (panel a), 40£o/' u O (panel b) and 70£o/^0 (panel c). The 
axes are in units of £o- The perturber is on the abscissa at position x = 0, and the flow enters the domain from the left. Only a limited 
region of the computational domain is shown (see section 4). The x (symmetry axis) and y (cylindrical radius) axes are labeled in units 
off . 
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Figure 5. Mass flux (in units of ^qPo^o, top), drag force (in units 
of £nPov%, centre) and gravitational force on the body (in units of 
^qPO^O' bottom) as a function of time (in units of £,o/vq). These 
parameters were computed from the results of the axisymmctric 
simulation (described in the text) using a spherical control volume 
of radius R = 5£o ■ 



the slope of the F d vs. R dependence appears to be higher 
in the numerical results than in the analytic model). 

In the bottom panel of Fig. 6 we show the gravitational 
force on the body due to the density structure of the nu- 
merical model. The gravitational force is larger than the net 
drag because momentum accretion onto the body produces 
an accelerating force. We also plot the gravitational force 
from the numerical simulation but excluding the contribu- 
tion from the dense wake behind the shock (see the lower 
sequence of points in the bottom frame of Fig. 6) . This force 
is comparable to the one obtained from the analytic density 
stratification given by Eq. (7), which only refers to the ma- 
terial in the free-streaming region before entry into the ax- 
ial wake. A comparison between the gravitational force with 
(upper sequence of points) and without (lower sequence) the 
material within the wake indicates that ~ 90% of the grav- 
itational force on the body comes from the material within 
the wake. 
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Figure 6. Time-averaged values of the net mass flux (in units 
of ^qPo^Oi t°p)i drag force (in units of JjQpo^o = (GM) 2 po/v 2 ,, 
centre) and gravitational force on the body (in the same units, 
bottom) as a function of the radius R of the control volume (R is 
given in units of £o). The dispersions arc indicated by the error 
bars. The gravitational force has been calculated in two ways: 
considering the contribution of all of the gas within the con- 
trol sphere (sequence of points on the upper part of the bottom 
frame), and eliminating the contribution of the material within 
the wake (lower sequence of points, bottom frame). The solid lines 
are the predictions from the analytic model described in section 
3. 



5 SUMMARY 

We present an analytic model for the flow generated by a 
point mass moving hypersonically within a homogeneous en- 
vironment. This model is based on the ballistic orbit theory 
(see, e. g., Bisnovatyi-Kogan et al. 1979), and is developed 
so as to obtain analytic expressions for the mass accretion 
rate and the non-linear gravitational drag. Since we include 
the contribution of the non-linear inner wake, there is no 
ambiguity in the definition of the minimum cut-off distance 
of the interaction, which turns out to be ~ 0.82^o. 

We find that the predicted mass accretion rate and grav- 
itational drag agree satisfactorily with the results from an 
axismmetric, isothermal simulation: 

• For the mass accretion rate we essentially find full agree- 
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ment between the analytic and numerical results (see Eq. 14 
and Fig. 6). This result in principle differs from the one of 
Moekel & Throop (2009), who obtain a significantly lower 
value for the accretion rate from their numerical simulation. 
The fact that we obtain a better agreement could be due to 
the considerably higher resolution of our simulation, or to 
the fact that we carry out a much longer time-integration 
(extending to ~ 70£o/vo, compared to ~ 1.5£o/«o for the 
simulation of Moekel & Throop 2009) . 

• For the net drag force Fd we obtain an agreement within 
~ 20 % between the prediction from the analytic model and 
the numerical simulations (see Fig. 6). Though the analytic 
and numerical drag forces have a reasonable quantitative 
agreement, it appears that the analytic model predicts a 
shallower Fd vs. R dependence (where R is the radius of 
the control volume enclosing the material assumed to pro- 
duce the drag) than the one obtained from the numerical 
simulation (see Fig. 6). 

There are several possible sources for this discrepancy be- 
tween the analytic and numerical drag forces. It appears 
that the limited numerical resolution of the simulation is 
not responsible for this effect, because we have repeated the 
simulation at 1/2 and 1/4 of the resolution (of the simula- 
tion presented in section 4) and obtain basically the same 
drag force. A possible source of the differences between the 
analytic and numerical Fd is the fact that the simulation 
has a finite Mach number (M = 5 for the upstream flow, 
see section 4), while the analytic model essentially has an 
infinite Mach number (i.e., zero gas pressure). Another dif- 
ference is that the numerical simulation has a rather broad 
wake region, while in the analytic solution it is assumed that 
the tail occupies a very narrow region surrounding the sym- 
metry axis. A third difference is that while in the analytic 
model the perturbed environmental gas effectively extends 
to infinity, the numerical simulation of course is carried out 
in a finite domain (see section 4). Given these clear differ- 
ences between the numerical and analytic models, the agree- 
ment that we find between the two can be regarded as quite 
successful. 

In this paper we have therefore derived an analytic 
recipe for the drag force F d (from a ballistic flow model), 
which is successfully reproduced by an axisymmetric numer- 
ical simulation. This recipe for Fd will be useful for carry- 
ing out simulations of compact bodies in motions influenced 
by gravitational drag. Possible examples are the motions 
of young stars within molecular clouds (see, e. g., Throop 
& Bally 2008; Chavarrfa et al. 2010), or the orbital decay 
of black holes in the centre of merging galaxies (Narayan 
2000; Escala et al. 2004, 2005; Dotti et al. 2006). Kim & 
Kim (2009) computed the nonlinear gravitational drag on a 
massive Plummer perturber in adiabatic axisymmetric sim- 
ulations and found that it is smaller than the linear theory 
predicts for supersonic bodies. This reduction of the drag 
force is accounted for correctly in our drag formula. 
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